A groupwise registration and tractography framework for cardiac myofiber architecture description by diffusion MRI: An application to the ventricular junctions

Background Knowledge of the normal myocardial–myocyte orientation could theoretically allow the definition of relevant quantitative biomarkers in clinical routine to diagnose heart pathologies. A whole heart diffusion tensor template representative of the global myofiber organization over species is therefore crucial for comparisons across populations. In this study, we developed a groupwise registration and tractography framework to resolve the global myofiber arrangement of large mammalian sheep hearts. To demonstrate the potential application of the proposed method, a novel description of sub-regions in the intraventricular septum is presented. Methods Three explanted sheep (ovine) hearts (size ~12×8×6 cm3, heart weight ~ 150 g) were perfused with contrast agent and fixative and imaged in a 9.4T magnet. A group-wise registration of high-resolution anatomical and diffusion-weighted images were performed to generate anatomical and diffusion tensor templates. Diffusion tensor metrics (eigenvalues, eigenvectors, fractional anisotropy …) were computed to provide a quantitative and spatially-resolved analysis of cardiac microstructure. Then tractography was performed using deterministic and probabilistic algorithms and used for different purposes: i) Visualization of myofiber architecture, ii) Segmentation of sub-area depicting the same fiber organization, iii) Seeding and Tract Editing. Finally, dissection was performed to confirm the existence of macroscopic structures identified in the diffusion tensor template. Results The template creation takes advantage of high-resolution anatomical and diffusion-weighted images obtained at an isotropic resolution of 150 μm and 600 μm respectively, covering ventricles and atria and providing information on the normal myocardial architecture. The diffusion metric distributions from the template were found close to the one of the individual samples validating the registration procedure. Small new sub-regions exhibiting spatially sharp variations in fiber orientation close to the junctions of the septum and ventricles were identified. Each substructure was defined and represented using streamlines. The existence of a fiber-bundles in the posterior junction was validated by anatomical dissection. A complex structural organization of the anterior junction in comparison to the posterior junction was evidenced by the high-resolution acquisition. Conclusions A new framework combining cardiac template generation and tractography was applied on the whole sheep heart. The framework can be used for anatomical investigation, characterization of microstructure and visualization of myofiber orientation across samples. Finally, a novel description of the ventricular junction in large mammalian sheep hearts was proposed.

ventricles and atria and providing information on the normal myocardial architecture. The diffusion metric distributions from the template were found close to the one of the individual samples validating the registration procedure. Small new sub-regions exhibiting spatially sharp variations in fiber orientation close to the junctions of the septum and ventricles were identified. Each substructure was defined and represented using streamlines. The existence of a fiber-bundles in the posterior junction was validated by anatomical dissection. A complex structural organization of the anterior junction in comparison to the posterior junction was evidenced by the high-resolution acquisition.

Conclusions
A new framework combining cardiac template generation and tractography was applied on the whole sheep heart. The framework can be used for anatomical investigation, characterization of microstructure and visualization of myofiber orientation across samples. Finally, a novel description of the ventricular junction in large mammalian sheep hearts was proposed.

Background
Mapping cardiac microstructure is essential to understand both the mechanical and electrical heart functions [1,2]. Although an active field of research, linking the organization of the cardiac microstructure to cardiovascular disease remains challenging. In particular, the mechanisms implicated in triggering and maintaining heart failure [3,4] or arrhythmias such as persistent AF [5] and ventricular fibrillation [6] are largely unknown.
Several cardiovascular imaging techniques such as micro-computed tomography (μCT) [7,8] and phase-contrast CT [9][10][11][12] combined with structure tensor image (STI) analysis have demonstrated their ability to map the cardiac fiber organization with a spatial resolution close to 10 μm isotropic [13]. They can provide information on myocardial-myocyte (myofiber) orientation and the myolaminar/sheetlet structure orientation. However, in vivo translation of phase-contrast X-ray imaging is limited by radiation exposure. Alternatively, destructive techniques like histology can now be performed in three-dimension at 1 μm isotropic in-plane resolution [14] and allow the assessment of cell-type distribution.
Cardiovascular diffusion tensor imaging (DTI) is a non-destructive method, allowing measurement of myocardial microstructure both ex vivo [15,16] and in vivo on humans [17] or animals models [18]. The clinical application of cardiac DTI was recently highlighted by several studies demonstrating how the variation of the microstructure between systole and diastole can be a marker of dilated cardiomyopathy [19] and hypertrophic cardiomyopathy [20]. Other in vivo studies [21,22] have also demonstrated that changes in diffusion tensor metrics, such as Fractional Anisotropy (FA), and Apparent Diffusion Coefficient (ADC) agree with late gadolinium enhancement (LGE) imaging and can depict areas with myocardial infarction and tissue remodeling. Despite growing progress in the field of MRI pulse sequences [23][24][25][26], hardware [27] and reconstruction methods [28], clinical in vivo DTI MRI remains limited in spatial resolution and signal-to-noise ratio (SNR) by clinical field strength, and scan times [29]. In vivo cardiovascular DTI must also take into account combined respiratory and cardiac motions, and compensated for them with a dedicated acquisition scheme [23]. In contrast, ex vivo cardiovascular DTI acquisition at ultra-high field on a static cardiac sample simplifies the acquisition scheme and allows for enhanced SNR, spatial and/or angular resolution. It therefore offers a better image quality compared to clinical acquisition providing an excellent research tool for anatomical investigation or validation of in vivo techniques.
Among all the ex-vivo studies using DTI, there is a large amount of evidence depicting the typical organization of the left ventricle (LV) myocardium characterized by the helix angle variation through the transmural wall [28][29][30][31]. Only a few studied describe sheetlets arrangement and fiber organization in both the LV and RV [30,31]. However, in our knowledge, there is a lack of a recent and high-resolution description of the microstructure organization of the junctions between the free wall of the left and right ventricles and the intraventricular septum (IVS) in large mammalian models often referred to as right ventricular insertion points (RVIP). A few decades ago, Kuribayashi and Robert [32] demonstrated the existence of abnormality of the myocardium at the junction associated with hypertrophic cardiomyopathy. By macroscopic inspection of tissue sections, the authors classified the histological characteristics of these junctions according to "myocardial disarray". Therefore, knowledge of normal fibers organization could be valuable for the identification of myocardial disarray and risk stratification [33].
The gold-standard method in image processing for creating a baseline map sharing spatial information with a group of subjects/samples (for instance a structural map describing the anatomy or the fiber organization and an associated parcellation) is named atlas or template. Such a framework normalizes the individual data to a common coordinate space and promotes the comparison between subjects of qualitative and quantitative information with one or multiple imaging modality sharing (or not) the same level of resolution. Spatial normalization is widely used in neuroimaging studies for structural and functional studies [34][35][36]. In cardiac microstructure studies, ex vivo DTI atlas has been initially proposed in 2006 on canine heart (0.3 × 0.3 × 0.9 mm 3 ) [37] and was applied later in humans (2× 2× 2 mm 3 ) [38,39], pigs [40], murine (0.043 × 0.043 × 0.043 mm 3 ) [41]. However, except for the murine study, these atlases have been created with low or anisotropic spatial resolution diffusion weighted (DW) images, and suffer from important voxel averaging of diffusion parameters.
In this work, we develop a template-based and tractography framework allowing us to resolve cardiac microstructural fiber anatomy of sheep hearts. An atlas is created with 3 ex vivo female sheep hearts fixed with high robustness sample preparation and acquisition protocol [16]. The template creation takes advantage of anatomical and diffusion-weighted images acquired at 9.4 T at an isotropic resolution of 150 μm, 600 μm respectively providing an accurate region-based analysis of myofiber architecture. We also emphasize the use of tractography to depict the global myofiber organization. To demonstrate the potential application of the proposed atlas and tractography visualization, we provide a novel description of sub-regions in the IVS and LV in which fibers present either an abrupt change in orientation or a notable deviation from the circular arrangement of the fibers depicted by the helix angle rule. The study focuses on: i) the presence of a singularity/fascicle of fibers at the junction of the posterior wall while noticing that the fiber arrangement mostly continues between both RV and LV wall and the septum. ii) the complex organization of the anterior junction where a division of the septal wall into the two fiber-bundles going through the aorta and pulmonary artery is visible iii) a singularity/fascicle close to the anterior junction going to the papillary muscle.

Sample preparation
The hearts of 4 female sheep (ovine) were explanted via sternal thoracotomy under general anesthesia, this protocol was approved by the Animal Research Ethics Committee (CEEA-050 Comité d'éthique pour l'expérimentation animale Bordeaux) in accordance with the European rules for animal experimentation. Three hearts were dedicated to MR experiments and one heart was dissected for macroscopic visualization. The three hearts (~12×8×6 cm 3 , heart weight = 150±10 g) were in relaxed state due to a first cardioplegic flushing before formalin fixation. They were then perfusion-fixed in 10% formaldehyde containing 2 mL of gadoterate meglumine, a gadolinium-based contrast agent to enhance the contrast between the myolaminae and cleavage planes and decrease the longitudinal relaxation time. To attenuate susceptibility artifacts at the border of the cardiac chamber, the hearts were removed from the solution and the cavity were filled with Fomblin, a perfuoropolyether with no 1 H signal when scanned using MRI. Lastly, whole hearts were immersed in Fomblin and sealed in a plastic container without rehydration. This procedure is described with additional details in Magat et al. [16].

Ex vivo anatomical and DWI acquisition
All MRI acquisitions have also been described previously in [16]. A summary is as follows: the experiments were performed at 9.4T (Bruker BioSpin MRI system, Ettlingen Germany) with an open bore access of 30 cm and a 200-mm inner diameter gradient (300 mT/m).
• A 3D FLASH sequence was applied to get anatomical images of the whole heart volume, at an isotropic resolution of 150 μm. The total acquisition time for each heart was 30 h 50 min. The resulting volume is referred to as the anatomical image. The parameters were: TE = 9 ms, TR = 30 ms, matrix size = 731 × 665 × 532, FOV = 110 × 100 × 80 mm 3 , flip angle = 30˚. An acceleration factor in the first phase dimension of 1.91 was used and the images were reconstructed with GRAPPA. The total acquisition time for each heart was 30 h 50 min.
• Diffusion-weighted (DW) images were acquired using a 3D diffusion-weighted spin-echo sequence at an isotropic resolution of 600 μm with 6 diffusion encoded directions, a single b-value of 1000 s/mm 2 and one b0 image. The parameters were: TE = 22 ms, TR = 500 ms, FOV = 100 × 80 × 110 mm 3 , matrix = 166 × 133 × 183. An accelerator factor in phase direction of 1.8 was used and the images were reconstructed with GRAPPA. The total acquisition time for each heart was 16 h 55 min.

Anatomical and DW images processing for template generation
The steps performed to generate the anatomical and the diffusion tensor template are detailed below and summarized in Fig 1. Unless otherwise mentioned, all registration steps and template generation were performed with ANTs library with the default parameters [42]. Estimation of diffusion tensor and derived tensor maps were performed using the MRtrix3 software [43,44]. All the processing was fully automatic, reproducible and executed using shell scripts without any input from the user with the exception of long axis alignment and ROIs and seeds definition. The corresponding transformations, ROIs and seed have been shared for reproducibility purpose (see data availability section).
A) Anatomical and FA template generation 1. Conversion from Bruker to Nifti format. Anatomical and DW images 2dseq files were converted in Nifti format with the dicomifier library (https://dicomifier.readthedocs.io/).
2. DW images pre-processing. Each of the DW images series was up-sampled by a factor of 2 using trilinear interpolation to reach a voxel size of 0.3×0.3×0.3 mm 3 . A denoising was performed using the MP-PCA method [45]. Finally, a B1 inhomogeneity field correction was applied with ITK-N4 [46]. Block diagram of the processing steps. A) Anatomical and DW images were first handled separately with a B1 inhomogeneity field correction. Then, anatomical images were manually aligned along the long axis using ITK-snap while diffusion tensor calculations were performed before any registration. Secondly, DW images were registered to the anatomical images. Model generation was performed using the symmetric normalization transformation (SyN) model using the anatomical and FA images. B) The individual diffusion tensors were then deformed into the template space and an average diffusion tensor was created. After this step, the diffusion 3. Diffusion tensor imaging estimation. Diffusion tensor calculations were performed prior any registration to avoid difficult transformation of the diffusion encoding matrix. Diffusion tensor maps (eigenvalues: λ 1 , λ 2 , λ 3 , ADC, FA, and color-coded FA (cFA) also known as Red-Green-Blue colormap) were obtained with MRtrix3. The first diffusion tensor eigenvector v1 corresponds to myofiber's main orientation [47,48]. The second v2 and third v3 eigenvectors were associated with sheetlet in-plane and normal directions, respectively.
4. Estimation of the cardiac mask. A binary mask was created and applied on both anatomical and DW images to segment myocardial tissue from regions with hyper intense signal corresponding to residual formaldehyde. Low and high cutoff thresholds were applied on the FA, trace, and DW images to define the binary mask as described previously in [16].
5. Co-registration of DW images on anatomical images. Anatomical and DW images were not always acquired during the same session because of the long scan time and scanner availability. Therefore, the alignment of each DW image on its respective anatomical images was needed and performed using Affine Registration using ANTs. A B1 homogeneity using ITK-N4 was also performed on anatomical images before the co-registration. Finally, FA maps were mapped to the space of the anatomical image as shown in Fig 1A. 6. Rigid Alignment of the anatomical and FA images. Images of each sample (anatomical and FA) were registered manually with ITK-snap [49] to align the long axis of the LV to the Z-axis of each volume and hence, obtain a similar alignment among hearts.
7. Template generation. The template generation [34], from coarse resolution to the final resolution, was performed iteratively with 4 levels of resolution (1 mm, 0.6 mm, 0.3 mm, 0.15 mm). Anatomical and FA images were used jointly using the following modality weights 1 × 0.5 respectively.
• At each level, the templates generated from the previous level were resampled to the current resolution and used as input to guide the current template creation using the standard script 'antsMultivariateTemplateConstruction2.sh' with the symmetric normalization transformation model (SyN) [50].
• For the first level (1 mm), the initialization sequence of template creation was similar to the one published in [51].
The same iterative strategy was performed to generate additional templates using only rigid transformation and affine transformation for comparison (Figs 2 and 3). After this step, the set of forward (or inverse) transformations for moving from individual subjects' space (also known as native space) to template space (or the inverse) were obtained. To determine registration quality, the Dice Similarity Coefficient [52] and Jaccard index [53] between each individual cardiac mask were computed for different transformation models (without transformation, after manual alignment, after Rigid, Affine and SyN registration). B) Diffusion tensor template generation 1. Mapping of individual diffusion tensor images to the template space. Individual diffusion tensor images were warped in the template space. This process includes both registration of the voxel at the new location and reorientation of the tensor via the 'antsRegistration' and 'ReorientTensorImage' command of ANTs. Preservation of eigenvector orientations related to the reorientation of the sample from the native space to the template space was carefully checked during all processing steps.
2. Diffusion tensor template generation and computation of derived maps. The transformed diffusion tensors were averaged via the AverageTensorImages' of ANTs to create the diffusion tensor template. Diffusion tensor maps were computed as described in the previous section A. 3.
C) Retrospective sampling to the standard in-vivo resolution.
1. The diffusion tensor images were retrospectively downsampled to the standard in-vivo resolution of 2.5×2.5×8mm 3 . The downsampling was performed in the log-Euclidean space using ANTs. Diffusion tensor maps were computed as described in the previous section A.3.

Tractography processing
Tractography was used for different purposes: i) Visualization of fiber organization, ii) Segmentation of sub-area depicting the same fiber organization, iii) Seeding and Tract Editing.

PLOS ONE
All the tractography steps described in this section were performed with MRtrix3 [43]. For each application, different algorithms were tested: • The deterministic tracking algorithm FACT [54] where streamline trajectory is determined by the principal eigenvector of the diffusion tensor.
Unless specified, we defined the tractography parameters as an FA-stopping threshold of 0.1, a maximum angle between steps of 60 degrees, a step size of 0.05 mm, a maximum length of 40 mm, a minimum length value of 1 mm. The choice of parameters, in particular, the termination criteria have no impact on the fiber orientation, and only affect the continuity of the streamlines. A brief reminder of the parameter impact is presented below: • The maximum angle between steps has a limited impact on connectivity as the propagation angle [57] (the angle between two adjacent vectors of the streamlines) is relatively low (<10˚) and therefore any value between 15˚and 60˚can be used.

PLOS ONE
has an impact on the streamline length: a higher value (0.2, 0.25) terminates streamlines in the lower anisotropy region.
• The criteria of maximum length also has an impact on the streamline generation. The streamline length does not reflect the size of the cardiomyocytes but it translates the overall arrangement of the myocytes.
Streamlines generated by FACT can be computed in the native space using the primary eigenvector of the individual diffusion tensors or in the template space using the primary eigenvector of the average diffusion tensor. The average diffusion tensor is depicting an averaged myofiber organization across samples. Streamlines generated by DET and PROB needed to be computed in the native space and warped in the template space to ease the comparison. Seeding ROI and tract filtering region of exclusion/inclusion were manually drawn in the template space and warped to the native space.
In the anterior and posterior junction, four cylindrical seeding regions (R1, R2, R3, R4) were manually drawn at the RVIP in the apex, mid-ventricular and base area. Based on the location of the region of interest (ROI) three pathways were arbitrarily defined to study the fiber organization at the junction: pathway A (from ROI R1 to R2) was the pathway inside the free wall of the RV. Pathway B (from ROI R1 to R2 to R3) described a pathway from the RV to the beginning of the free wall of the LV while pathway C (from ROI R1 to R2 to R4) described a pathway from the RV to the IVS. For each pathway, streamlines were computed and filtered by identifying tracts that pass through all the corresponding ROIs. For each tracking algorithm tested, the resulting number of streamlines were counted.
To investigate the distribution of propagation angle (PA), a tractography-based metrics quantifying the curvature of streamlines, a new tractogram was generated by FACT for samples #1, #2, #3 in the posterior RVIP. We defined the tractography parameters as an FA-stopping threshold varying from 0.1 to 0.25, a maximum angle between steps of 45 degrees, a step size of 0.3 mm, a maximum length varying from 10 mm to 40 mm, and a minimum length value of 1 mm. The propagation angle was computed along myofiber trajectories within the principal eigenvector field similarly to Mekkaoui [56] and mapped as a scalar volume.

Dissection study and macroscopic anatomy visualization
A heart dedicated to macroscopic examination was added to the study to confirm the presence of this substructure. The dissection was performed with knowledge of the identified structure and 3D rendering visualization was used to guide the dissection. Macroscopic examination was performed as follows: • The LV, RV and then the IVS were transected from base to apex to remove the anterior part of the heart and isolate the posterior wall of the heart.
• Transverse (short-axis) cuts of the sample were performed progressively from base-ventricular level to mid-ventricular level.
• A coronal cut (long axis) was performed on line ranging from the RV cavities to the LV posterior wall going through the suspected region.
Photos were taken with a Canon PowerShot SX710 HS and a Samsung Galaxy SM-G960F with the surgical light of the operating room. Rigid registration was manually performed using ITK-snap to align photography on anatomical images.

Results
Fig 2 shows sagittal, coronal and transverse views of the anatomical and FA templates with different transformation models used during the template creation. The template includes not only the ventricles but also the basal parts of the heart, including the valve, atria, the aorta or the pulmonary artery. Misregistration areas were more visible with rigid or affine transformations than with the SyN transformation model. Several differences are evident through visual inspection: (i) a sharper definition of the myocardial wall (Figs 2 and 3A-3F); (ii) an enhanced definition of the division of the His bundle in left and right bundles ( Fig 3B); (iii) a more accurate delineation of the leaflets of the tricuspid valve (Fig 3B, 3E and 3F). Consistently misregistered areas using the SyN transformation model were noticeable in the basal part (Fig 2) and for most of the small or large vessels ( Fig 3E) going through the myocardium or around (like the coronary sinus) due to the variability of the vessel locations. S1A and S1B  #1, #2, #3, on the average diffusion tensor images, and all together, respectively). All datasets were warped in a template space to ease the comparison. In the short-axis view, the blue color indicates that fiber orientation in the papillary muscle or in the IVS is parallel to the long axis on the heart. On the other hand, green and red colors indicate the orientation of the fibers in the short axis plane in Anterior-Posterior (AP) or Left-Right (LR) orientation respectively. While a smooth gradient of color from green to red (or vice versa) indicates a normal arrangement of the cardiac myofibers, the existence of small sub-regions in the LV indicates either a rapid change in orientation or apparent boundaries. Especially: fiber arrangements in Inferior-Superior (IS) orientation are visible: i) from the basal area to the mid-ventricular level close to the papillary muscle, ii) at the intersection of the LV and RV from the basal area to mid-ventricular level, iii) in the basal area of the IVS, divided into two fiber-bundles that go, on the one hand, to the pulmonary artery and, on the other hand, to the aorta.
To further investigate the existence of abrupt changes in orientation, we first focused our analysis on the posterior wall. Fig 4A-4C displays coronal, sagittal and transverse views of the anatomical and FA templates at the junction of the RV, LV and IVS in the basal area. Fiber orientation can be visually identified on the anatomical template images due to 150 μm resolution and the T2 � contrast. A brighter image intensity is noticeable in the layer of the IVS close to the RV and at the junction of the ventricles. The FA derived from the diffusion tensor template is relatively homogeneous except for two regions with lower values. They are i) a line dividing the IVS into two areas and ii) a triangle located at the junction. Fig 4D and 4E displays the cFA Anatomical images of the IVS and the picture of the dissected sample are shown in transversal orientation in Fig 6A and 6B, respectively. An overlay of both images is shown in Fig 6C after rigid registration and shows an acceptable agreement by visual inspection. The coronal cut in the posterior area is indicated by a yellow line ranging from the RV cavities to the LV in Fig 6D. A section of the dissected sample is presented using five pictures for which the sample has been slightly rotated (Fig 6E-6I). Images reveal macroscopic features depicting the fiber orientation changes and support the evidence of edges. As a comparison, streamlines are displayed in short axis and long axis orientation in Fig 6J and 6K respectively. The fiber orientation display between Fig 6K and 6G or 6L shows similar fiber orientation. Fig 7A and 7B displays the cFA maps of the anterior junction close to the base in axial and sagittal views. The white dash line indicates the plane of the sagittal view. A smooth transition of fiber orientation is visible in the LV and the endocardium part of the IVS (with a gradient of green to purple). In the epicardium part of the IVS, we observed a division of the fibers inside the IVS into the fiber-bundles going through the aorta (yellow arrow) and the pulmonary artery (purple arrow). The third fascicle of fibers (in pink) located close to the endocardium in the LV was also noticeable with an orientation differing from the surrounding structures. For each mentioned area, streamlines were computed using seeds on the averaged diffusion tensor with the FACT algorithm. To facilitate the rendering, the sample is flipped and the wall of the pulmonary artery and the aorta is manually segmented and plotted as a surface mesh in Fig  7C-7H. Tractography shows a continuity between the IVS and the wall of the pulmonary artery (in green Fig 7F and 7I), a continuity between the IVS and the wall of the aorta (in green Fig 7G and 7I) and the existence of fibers that split into two branches: one going to the papillary muscle, the other surrounding the top part of the LV (in red Fig 7H and 7I). Some abruptly between adjacent voxels in the IVS and in the middle of the RVIP (yellow arrow) depicting a triangular shape in the coronal view. (E) streamlines with the color-coding of the cFA maps (D). Note the presence of a channel or fascicles of streamlines going from the RV to the IVS (in green). https://doi.org/10.1371/journal.pone.0271279.g004

Fig 5. Streamlines of the posterior singularity with different tracking techniques on template and individual diffusion tensors.
For each representation, the anatomical template is inserted for reference. Top) Streamlines were computed using a masking region of interest and the FACT algorithm on the primary eigenvector of the average tensor in template space and displayed in coronal (A) and sagittal view (B). Down) A different tracking technique is shown. Few seeds were manually defined, then streamlines were generated for sample #1 (C), #2 (D), #3 (E) using the Tensor_Det algorithm on each individual tensor in native space and warped into template space without ROI restriction. streamlines going from PA to the LV and from AO to the RV (yellow arrows) are also found close to the epicardium of the ventricles.
To illustrate the added value of the seed-based tractography and analyze the global myofiber architecture, three arbitrary pathways of fibers are defined in the anterior and posterior junctions in the mid-ventricular area. The resulting streamlines for sample #1 are proposed in Fig  8 in 3D rendering views. In the posterior wall (down), the connection of fibers originating from the RV into the LV free wall and the IVS is detectable (note that the unidirectional direction is arbitrary) (pathway A). Existing pathways (B) and (C) were found connecting. In the anterior wall (top), streamlines split into two output tracks but the fascicle of fibers in the IVS leaves the image plane and goes to the wall of the pulmonary artery as described in

PLOS ONE
free wall (pathway A). Existing pathways (B) and (C) are found connecting. On the contrary, in the anterior wall (top), the connection of fiber originating from the RV into the LV free wall is limited to output tracks located on the epicardium surface of the sample and the pathway (C) is not found. Table 1 summarizes the streamline frequency percentages for the 3 samples in apex, mid-ventricular, base areas for the three tracking algorithms. The ratios are all equal to 100% for pathway A and close to 100% for pathway B except for the posterior basal area (69.2%, 2.1%, 62.3% using the DET algorithm for samples #1, #2, #3 respectively) and the anterior apex area (100%, 32.1%, 2.6% using the DET algorithm for samples #1, #2,#3 respectively). For pathway C, the ratios are close to 0% in the anterior junction (blue square) and range from 0% to 100% in the posterior junction (red square). The DET and PROB algorithms gave approximately the same ratio, the use of the FACT algorithm tends to decrease the ratio by a factor of 2 or more. Variability over samples is noticeable for pathway B in the posterior mid area and the anterior apex area and for pathway C in the posterior basal and apex area. The corresponding ratios are 100%, 2.4%, 100% and 13.6% 0% 97.1% for samples #1, #2, #3. Table 1. Global myofiber architecture of the junction using seed-based tractography. Streamlines count as a percentage of the 3 samples in apex, mid-ventricular, base areas for the 3 different tracking algorithms. The ratio is defined as the number of resulting streamlines over the total number of requested streamlines. Variation in connectivity and variability over samples are particularly noticeable for the pathway C. Nevertheless, in all cases, the ratios were found inferior for pathway C for all areas but remained close to 100% (square with red dotted line) in the posterior junction for the base and mid-ventricular level while values were close to 0% in the anterior junction (square with blue dotted line).

Discussion
Using anatomical and DW images, we proposed a new framework combining atlas creation and tractography to characterize the cardiac microstructure of ex-vivo sheep hearts. We first built a 150 μm isotropic anatomical template and a 300 μm isotropic voxel size FA template. Assessment of the accuracy of inter-subject spatial normalization was done qualitatively by visual inspection. Although misregistration of vessels was visible in different areas of the myocardium, the proposed pyramidal approach to template creation resulted in excellent anatomical consistency.
The main novelties of the proposed template were to cover both ventricles, the atria, and a fraction of the arteries at an unprecedented spatial resolution of 600 μm isotropic for the DTI acquisition. The diffusion metric distributions were found close to the individual sample (S2 Fig) and agree with prior studies on other species [37,41] ensuring that the proposed DTI template creation framework preserves the quantitative DTI metrics. Although image-based quantitative measurements were made to assess the accuracy of the model (S1 Fig), the preservation of the FA distribution is also a strong argument. In case of mis-registration, tensor averaging would have decreased the anisotropy [58].
As shown in Fig 4, spatial variation of the FA was observed in individuals and averaged diffusion tensors. This could indicate i) intra-voxel variability with an abrupt change of fiber orientation within the voxel leading to reduced anisotropy (in the mid-wall of the septum) ii) the presence of inter-subject anatomy variability (close to the edge of the posterior singularity in the averaged diffusion tensor) iii) the potential existence of registration errors at the voxel scale. A combination of all three scenarios is also likely probable.
The study emphasizes the use of tractography to depict the cardiac fiber organization. This approach is commonly used in neuroscience to study brain connectivity [59][60][61]. Except for the studies mentioned by Sosnovik et al. [62], its use in cardiac diffusion remains marginal and focuses mostly on whole heart tractography for visualization purposes. It is important to highlight that there is an order of magnitude between the size of cardiomyocytes (around 120 μm × 40 μm in diameter [63]) and both the minimum length of the computed streamline (1 mm) and also the voxel size of the reconstructed image (300 × 300 × 300 μm 3 ). In this work, the streamlines were used to give an insight of the myofiber architecture, Figs 5-7 show how this 3D representation depicts both the local variation of the fiber orientation but also the continuity of the fiber in the surrounding area. The continuity is dependent on the DW analysis methods, the tracking algorithms and the input parameters [64]. In this study, a comparison of diffusion-weighting techniques was not possible due to the unique b-value and the acquisition of 6 diffusion directions. We compared only the use of three different tracking algorithms. The angle parameter of 60˚was found to be a maximum to maintain the physical behavior of the fiber (an angle of 90˚or more results in undulating streamlines). The cut-off parameter was set to 0.1 to maximize connectivity as shown in S4 Fig. The max length parameter was easier to set either depending on the expected size of the object or the desired rendering. Pitfalls of tractography techniques are numerous [65] and have been raised in many studies even with state-ofthe-art algorithms [59,66]. Nevertheless, investigation of the spatial continuity was found difficult using only a visual assessment of the eigenvectors on 2D slices or with a Multiplanar reformation or reconstruction (MPR) viewer.
In our work, the first limitation is the low angular resolution of the DW dataset. A larger number of diffusion directions is recommended for a better estimation of the diffusion tensor [67]. Recent ex vivo cardiac studies reported a number of gradient directions from 7 [68], 30 [69], 32 [70] to 61 [68]. Nevertheless, as demonstrated in [68], the trade-off between signal-tonoise ratio and angular resolution are not the only parameters to take into account. The heterogeneity of diffusion directions is better depicted using a high spatial resolution and without voxel anisotropy limiting the voxel averaging of the DTI metrics. By focusing on the spatial resolution and SNR, the proposed approach excludes the acquisition of high angular sampling of the Q-space in regard to the long scan time and prevents from applying more advanced diffusion characterization models like q-ball imaging (QBI) or diffusion spectrum magnetic resonance imaging (DSI).
The second limitation is the tensor model that constrains the fiber orientation to one unique direction per voxel. DSI tractography could provide a more realistic representation of the myofiber architecture [62]. The democratization of undersampled acquisition [28] in preclinical settings could also be a potential solution to increase the angular resolution while limiting the scan time.
To demonstrate the added value of the tractography, the structural organization of the anterior and posterior junction was characterized. The study focuses on one fiber-bundle in the posterior junction and three fiber-bundles in the anterior junction. Tractography in a selected ROI was useful to visualize the spatial coverage of the region independently of the surrounding environment ( Fig 5A). Tractography with a specific region seeding (e.g., the user manually delineates only a small ROI in the fiber bundle) was useful to display the fiber organization without apriori knowledge of the global myofiber arrangement (Figs 5C-5E and 7E-7H). In the last part of the study (Figs 8 and S5), generated streamlines were constrained by a user rule and must traverse all specified ROIs to be included in the final tractogram. The ROI was defined on both sides of the junction at three different levels of the ventricle (apex, mid, base). Note that the tested rule did allow crossing pathways between levels and represent a fraction of the structural organization. We then analyzed qualitatively and quantitatively the number of generated streamlines that respect a user rule. A ratio of 100% indicating that all requested streamlines have been generated by the tracking algorithms. Therefore the resulting value is not only dependent on the tracking parameters but also on the request number of streamlines. In all cases, the ratios were found inferior for pathway C for all areas but remained close to 100% in the posterior junction for the base and mid-ventricular level while values were close to 0% in the anterior junction indicating the absence of connectivity. We also analyze tractography-based metrics to assess the curvature of streamline also referred as the propagation angle. As demonstrated by Mekkaoui [56], the existence of a high propagation angle could be correlated with changes in electrical activation waves. The propagation angle displayed in S7 Fig is coherent with the fiber orientation discontinuity in the junction. A homogenous low value is mostly found in the myocardium except where fiber discontinuity occurs (in the IVS and at the RVIP) where it is slightly increased. The presence of abrupt change of fiber orientation may be associated with a conduction block or ectopic beat in pathological hearts. Combination of submillimeter diffusion weighted acquisition with dedicated high-level image processing could increase our knowledge of cardiac microstructure and help to understand cardiac electrophysiology by differentiating normal to remodeled or pathologic discontinuities [71].
The spatial resolution of 0.6 mm isotropic, and the acceptable inter-variability over samples, enabled the identification of fibers in the IVS and LV that present either a rapid change in orientation or some deviation or a notable deviation from the circular arrangement of the fiber depicted by the helix angle rule. At first glance, we noticed a more complex structural organization of the anterior junction in comparison to the posterior junction. The global myofiber organization was found closely linked to the surrounding structure including the pulmonary and the aorta. The "Y" pattern, depicted in Fig 7B, is also associated with the location of the anterior-septal perforator arteries recently mentioned in [72] referred as ventricular myocardial extensions [73]. Finally, visual inspection and streamline connectivity of the anterior junction suggest that most fibers end in the anterior wall while a strong connectivity fiber is visible in the posterior wall. Indeed, a smoother transition of fiber orientation was noticed in the posterior wall (gradient of green to purple in Fig 4) in the endocardial part of the LV while fiber orientation changes abruptly between adjacent voxels in the IVS in agreement with [74] and the middle of the RVIP (blue in Fig 4) depicting a triangular shape in the coronal view. The macroscopic size (1 cm 3 volume) of this substructure suggests that experimental observation of the myofiber orientation could be seen by the eye without magnification. We proceeded to a dissection of the heart of another animal that confirmed the presence of this substructure. However, the macroscopic examination is a destructive technique whose cut might alter the integrity of the revealed structure. This issue is quite common in histology during tissue preparation with fixation and slicing. It can lead to misalignment [75], shrinkage, distortion [76] however excellent correlation between standard histology and DTI have been described in the literature for cardiac fiber orientation [76].
While clinical implication is out-of-scope of this study, it is interesting to notice that uncommon idiopathic ventricular arrhythmias have been reported in the both anterior RVIP region [77,78] and posterior RVIP region [79][80][81] also mentionned the basal inferoseptal left ventricle and could be related to the specific myofiber organization of these region. Additionally, Domenech-Ximenos et al. [82] have also recently reported that highly trained endurance athletes present a higher prevalence of focal LGE with a triangular pattern in particular at the RV attaches to the septum, a visual similarity between the identified substructure defined as posterior singularity and the area defined by contrast enhancement is noticeable (Figs 1 and  2). Different hypothesis could be envisioned: i) the presence of boundaries due to the drastic change in fiber orientation might be more sensitive to structural remodeling. ii) the top/down orientation of the fiber could be more sensitive to contraction-induced mechanical stress. Nevertheless, the prevalence of the spotty pattern of LGE localized at the junction of the right ventricle to the septum is not associated with life-threatening arrhythmias and sudden cardiac death in the athlete [83]. To investigate the possibility of in-vivo imaging of the identified structure, the diffusion tensor of samples #1, #2¸#3, and the average diffusion tensor were downsampled to 2.5×2.5×8mm 3 . The new voxel size (50 mm 3 ) is 231 times larger than the previous one (0.216 mm 3 ) and increases the partial volume effect. Nevertheless, a few remaining voxels (<10 per sample) displayed in S6 Fig still shows characteristic fiber orientation in the RVIP. However, shifting to in-vivo acquisition while maintaining this image quality will remain extremely challenging due to the combined cardiac and respiratory motion.
The finding must be balanced for two different reasons: i) the fix state of the ex-vivo samples and ii) the low number of datasets. i) Hearts cardiac muscle is fully relaxed due to cardioplegic flushing before formalin fixation. Therefore, hearts were imaged in a relaxed state but were not in a complete diastolic phase. A direct comparison with in vivo DTI measures is complicated because the state of ex vivo hearts does not match the in vivo contractile state [16,48]. Therefore, the finding cannot be directly transposed to in vivo where the fiber angles vary throughout the cardiac cycle. However, limited changes are observed for fiber organization [18,20,84] between in vivo cardiac phase (systolic and diastolic) and ex vivo states (relaxed and contracted) and only sheetlets orientations are highly dependent of contractile state. Thus, this fiber-bundle which is aligned in the ex vivo state perpendicular to the usual fiber orientation is probably present in vivo at any contraction stage. ii) The use of three hearts only is a limitation since it does not allows us to explore the inter-subject variability. However, our study focuses on a macroscopic area which present in the three individuals. Moreover, efforts were dedicated to limit the inter subject variability by using animals of identical strain, age and sex as well as by using the exact same sample preparation protocol. As a comparison the first atlas of the human heart was built with 10 hearts [37] while using individuals of different ages and sex which could induce important inter-subject variability. That is the reason why this number appear sufficient to enforce the conclusion of our study and thus was kept low for ethical consideration.
Finally, mechanical and functional aspects are not discussed in the current study but the individual DW images, diffusion tensors, averaged diffusion tensors of the junction are already available in the Zenodo platform and delivered with shell script for computing the diffusion tensor maps and the streamlines. The conversion of the eigenvectors to CSV format for mesh generation and integration into EP/mechanical models is also possible (see Data availability section) and the visualization of the results can be done with well-documented open-source software.

Conclusion
Demonstrating that the myocardial-myocyte orientation or myolaminar/sheetlet structure changes are correlated with pathology could theoretically allow the definition of relevant quantitative markers in clinical routine to diagnose heart pathologies. A whole heart diffusion tensor template representative of the global myofiber organization over species is therefore crucial for comparisons across populations. In this paper, we proposed a new framework combining template creation and tractography to characterize the cardiac microstructure of three ex vivo sheep hearts and demonstrate the potential application by providing a novel description of the ventricular junction in large mammalian sheep.